home *** CD-ROM | disk | FTP | other *** search
/ The Fatted Calf / The Fatted Calf.iso / Applications / Graphics / MapMaker / Source / proj.c < prev    next >
Encoding:
Text File  |  1990-12-04  |  9.0 KB  |  392 lines
  1. /*
  2. ** proj.c
  3. ** map projection calculation and generation routines
  4. ** for project : MapMaker  
  5. ** using NeXTStep and mach Unix.
  6. ** CPSC 414 Assignment No. 4 Project
  7. **  
  8. ** Programmed by Bradley Head and Thomas Burkholder 
  9. ** December 1990
  10. */
  11.  
  12. #import "proj.h"
  13. #import "pointdata.h"
  14. #import "appkit/graphics.h"
  15.  
  16. /********************* THE MAP PROJECTIONS (Forward Sol'n) **********************************/
  17.  
  18. float correct(float f)
  19. {
  20.     if (f>=PI)
  21.         return f-(2*PI);
  22.     if (f < (-(PI)))
  23.         return f+(2*PI);
  24.     return f;
  25. }
  26.  
  27. /*
  28. **  Map Projection conversion functions 
  29. */
  30.  
  31. int Cylindrical(Point *posn, ProjParam *init, Point *pp)
  32. { 
  33. float coef;
  34.  
  35. coef = init->radius*cos(init->phi1);
  36. pp->x = (posn->x - init->lon0)*coef;
  37. pp->y = init->radius*sin(posn->y)/cos(init->phi1);
  38. if (pp->x > (PI*coef)) {
  39.     pp->x -= ((2*PI)*coef);
  40.     return 1;
  41.     }
  42. if (pp->x < ((-(PI))*coef)) {
  43.     pp->x += ((2*PI)*coef);
  44.     return 1;
  45.     }
  46. return 0;
  47. }
  48.  
  49.  
  50. int EckertIV(Point *posn, ProjParam *init, Point *pp)
  51. { 
  52. float coef;
  53.  
  54. coef = .42224*init->radius*(1 + cos(posn->y - init->phi1));
  55. pp->x = coef*(posn->x - init->lon0);
  56. pp->y = 1.32650*init->radius*sin(posn->y - init->phi1)/cos(init->phi1);
  57.  
  58. if (pp->x  > (PI*coef)) {
  59.     pp->x -= ((2*PI)*coef);
  60.     return 1;
  61.     }
  62. if (pp->x < ((-(PI))*coef)) {
  63.     pp->x += ((2*PI)*coef);
  64.     return 1;
  65.     } 
  66. return 0;
  67. }
  68.  
  69. int Mercator(Point *posn, ProjParam *init, Point *pp)
  70. {
  71. float coef;
  72.  
  73. coef = init->radius*0.6;            /* 0.6 scales to fit into output window with radius = 1.0 */
  74. pp->x =coef*(posn->x - init->lon0);
  75. pp->y = init->radius*log(tan(.785398 + posn->y/2));
  76. if (pp->x > (PI*coef)) {
  77.     pp->x -= ((2*PI)*coef);
  78.     return 1;
  79.     }
  80. if (pp->x < ((-(PI))*coef)) {
  81.     pp->x += ((2*PI)*coef);
  82.     return 1;
  83.     }
  84. return 0;
  85. }
  86.  
  87. int Mollweide(Point *posn, ProjParam *init, Point *pp)
  88. { 
  89. float coef;
  90.  
  91. coef = .9003*init->radius*cos(posn->y - init->phi1);
  92. pp->x = (posn->x - init->lon0)*coef;
  93. pp->y = 1.4142*init->radius*sin(posn->y - init->phi1);
  94. if (pp->x > (PI*coef)) {
  95.     pp->x -= ((2*PI)*coef);
  96.     return 1;
  97.     }
  98. if (pp->x < ((-(PI))*coef)) {
  99.     pp->x += ((2*PI)*coef);
  100.     return 1;
  101.     }
  102. return 0;
  103. } 
  104.  
  105. int Ortho(Point *posn, ProjParam *init, Point *pp)
  106. {
  107. float cos_c;
  108. float dl; 
  109. dl = posn->x  - init->lon0;
  110.  
  111. pp->x = init->radius*cos(posn->y)*sin(dl); 
  112. if ((posn->y != -PI) && (posn->y != PI))
  113.     pp->y = init->radius*(cos(init->phi1)*sin(posn->y) - 
  114.     sin(init->phi1)*cos(posn->y)*cos(dl));
  115. else
  116.     pp->y = init->radius*cos(posn->y)*cos(dl);
  117. if (posn->y == PI)                            /* conditionals for line removal */
  118.     pp->y = - pp->y;
  119. cos_c = sin(init->phi1)*sin(posn->y)+cos(init->phi1)*cos(posn->y)*cos(dl);
  120. if (cos_c < 0)
  121.     pp->pen = SKIP;
  122. return 0;
  123. } 
  124.  
  125.  
  126. int Sinusoidal(Point *posn, ProjParam *init, Point *pp)
  127. {
  128. float coef;
  129. float xofs;
  130.  
  131. coef = init->radius*cos(posn->y - init->phi1);
  132. pp->x = (posn->x - init->lon0)*coef ; 
  133. pp->y = init->radius*(posn->y - init->phi1);
  134. if (pp->x > (PI*coef)) {            /* clipping conditonals for animation */
  135.     pp->x -= ((2*PI)*coef);
  136.     return 1;
  137.     }
  138. if (pp->x < ((-(PI))*coef)) {        /* conditionals for animation */
  139.     pp->x += ((2*PI)*coef);
  140.     return 1;
  141.     }
  142. return 0;
  143. }
  144.  
  145.  
  146. int Stereo(Point *posn, ProjParam *init, Point *pp)
  147. { 
  148. float dl,k,r; 
  149. float Len;
  150. dl = posn->x - init->lon0;
  151. r = init->radius*0.6;
  152. k = 2/(1+sin(init->phi1)*sin(posn->y)+cos(init->phi1)*cos(posn->y)*cos(dl));
  153. pp->x = r*k*cos(posn->y)*sin(dl);
  154. pp->y = r*k*(cos(init->phi1)*sin(posn->y) - sin(init->phi1)*cos(posn->y)*cos(dl));
  155. Len = sqrt((pp->x)*(pp->x) + (pp->y)*(pp->y));
  156. if (Len > init->radius*1.2)
  157.     pp->pen = SKIP;
  158. return 0;
  159. }
  160.  
  161. int None(Point *posn, ProjParam *init, Point *pp)
  162. { 
  163. float dl; 
  164. float Len;
  165. float coef;
  166. coef = init->radius;
  167. dl = posn->x - init->lon0;
  168. pp->x = coef*dl;
  169. pp->y =  coef*(posn->y - init->phi1);
  170. if (pp->x > (PI*coef)) {            /* clipping condtionals for animation */
  171.     pp->x -= ((2*PI)*coef);
  172.     return 1;
  173.     }
  174. if (pp->x < ((-(PI))*coef)) {
  175.     pp->x += ((2*PI)*coef);
  176.     return 1;
  177.     }
  178. return 0;
  179. }
  180.  
  181.  /****************************************************/ 
  182.  
  183.  void DoFrame(PointList *framelist) {
  184.  /*
  185.  **  This function creates two frame lines 
  186.  **   for the projections
  187.  **   all except stereographic and orthographic
  188.  */        
  189.   float lat, lon;
  190.  Point framept;
  191.  framept.pencolour = NX_DKGRAY;
  192.    for (lon = -180;lon <= 180;lon += 360) {
  193.      for (lat = -75;lat <= 75;lat += 5)  {
  194.         if (lat == 75)
  195.             framept.pen = UP;
  196.         else
  197.             framept.pen = DOWN;
  198.         framept.y = lat *toRADS;
  199.         framept.x = lon *toRADS;
  200.         addToPointList(framelist ,&framept);
  201.         }
  202.     }
  203.  }
  204.  
  205.  void  drawCircleFrame(float radius, PointList *pout) {
  206.  /* 
  207.  ** This function generates
  208.  ** a framing circle
  209.  ** around stereographic and
  210.  ** the orthographic
  211.  ** projections.
  212.  */
  213.  Point p;
  214.  float  step, x;
  215.  p.pencolour = NX_DKGRAY;
  216.  step = radius/100.0;
  217.  for (x = -radius;x <= radius;x += step) {    /* to top half of circle */
  218.      p.y = sqrt(radius*radius - x*x);
  219.     p.x = x;
  220.     p.pen = DOWN;
  221.     addToPointList(pout,&p);
  222.     }
  223.  for (x = radius-step;x >= (-radius+step); x -= step) {  /* do bottom half */
  224.     p.y = -(sqrt(radius*radius - x*x));
  225.     p.x = x;
  226.     p.pen = DOWN;
  227.     addToPointList(pout,&p);
  228.     }
  229. p.y = 0.0;
  230. p.x = -radius;
  231. p.pen = UP;
  232. addToPointList(pout,&p); 
  233.  }
  234.  
  235. /****************************************************/
  236.  
  237. void DoGrid(PointList *gridlist) {
  238. /*
  239. ** This function generates the graticules
  240. ** for the output view 
  241. ** for all the projections
  242. ** it increments at 15 degree intervals
  243. */
  244. int lat, lon;
  245. Point  gridpt ;
  246.  
  247. gridpt.pencolour = NX_LTGRAY;
  248. gridpt.pen = DOWN;
  249. /*  
  250. ** Compute the lines of latitude
  251. */
  252.  
  253. for (lat = -75;lat <= 75;lat += 15) { 
  254.     for (lon = -180;lon <= 180; lon += 5)
  255.         { 
  256.         if (lon == 180 )
  257.             gridpt.pen = UP;
  258.         else 
  259.             gridpt.pen = DOWN;
  260.         gridpt.y = lat *toRADS;
  261.         gridpt.x = lon *toRADS;
  262.         addToPointList(gridlist,&gridpt);
  263.         }
  264.     } 
  265. /*  
  266. ** Now compute the Meridians
  267. */
  268.  
  269. for (lon =  -180;lon < 180; lon += 15) {
  270.     for (lat = -75;lat <= 75;lat += 5 ) 
  271.         {
  272.         if (lat == 75)
  273.             gridpt.pen = UP;
  274.         else
  275.             gridpt.pen = DOWN;
  276.         gridpt.y = lat *toRADS;
  277.         gridpt.x = lon *toRADS;
  278.         addToPointList(gridlist,&gridpt);
  279.         }
  280.     } 
  281. }
  282.  
  283. void ComputePoints(int (*proj)(Point *posn, ProjParam *init, Point *pp), ProjParam *init,PointList *inlist,PointList *outlist)
  284. {
  285. /* 
  286. ** This function computes the output
  287. ** points given the input list of  (lat,lon) pairs
  288. ** it  calls the appropriate map projection routine
  289. ** determines line clipping in the form of skipping
  290. ** unnecessary points (points that won't be drawn to output)
  291. ** It places the computed points in an outlist list of points
  292. ** that can then be plotted
  293. */
  294. int i ;
  295. int oldsplit,split;
  296. Point pp, old;
  297. Point *pt;
  298. int a =1;
  299. old.pen = UP;
  300. oldsplit = 1;
  301. for (i = gotoPointInList(inlist,0,&pt);(i);i=gotoNextPointInList(inlist,&pt))    
  302.     {
  303.     
  304.     pp = *pt;
  305.     split = (*proj)(pt,init,&pp);
  306.     if ((old.pen == DOWN) && (pp.pen == SKIP)) 
  307.         old.pen = UP;
  308.     if (split  != oldsplit)
  309.         old.pen = UP;
  310.     if (a) 
  311.         a = 0;
  312.     else
  313.         addToPointList(outlist,&old);
  314.     old = pp;
  315.     oldsplit = split;
  316.     }
  317. if (inlist->quantity) {
  318.     pp.pen = UP;
  319.     addToPointList(outlist,&pp);
  320.     }
  321. }
  322.  
  323. int convertPoints(PointList *pin, PointList *pout, int gridon, ProjParam *init)
  324. {
  325. /* This function is the externally called function
  326. ** It handles whether to compute and draw the 
  327. ** grid and frames and computes the output points
  328. ** by calliing ComputePoints with the appropriate 
  329. ** map projection.
  330. */
  331. PointList grid, frame;        /* point list for the frame around the projection */
  332. ProjParam frameparam;        /* frame drawing initialization paramaters */
  333.  
  334. frameparam.radius  = init->radius;    /* get the current radius */
  335. frameparam.lon0 = 0.0;            /* fix the central meridian for frame computation */
  336. frameparam.phi1 = init->phi1;    /* get the current central line of latitude */
  337.  
  338. init->lon0 = correct(init->lon0);    /* correct the offset for animation */
  339. newPointList(&grid);            /* generate a new list for the grid */
  340. newPointList(&frame);            /* generate a new list for the frame */
  341. DoFrame(&frame);                /* compute the frame */
  342. if (gridon)                        
  343.     DoGrid(&grid);            /* if want the grid displayed then compute the grid */
  344.  
  345. switch (init->proj)                /* select the appropriate projection to compute */
  346.     {
  347.     case CYLINDER: 
  348.         ComputePoints(&Cylindrical,init,&grid,pout);
  349.         ComputePoints(&Cylindrical,&frameparam,&frame,pout);
  350.         ComputePoints(&Cylindrical,init,pin,pout);
  351.         break;
  352.     case ECKERT:  
  353.         ComputePoints(&EckertIV,init,&grid,pout);
  354.         ComputePoints(&EckertIV,&frameparam,&frame,pout);
  355.         ComputePoints(&EckertIV,init,pin,pout);
  356.         break;
  357.     case MERCATOR: 
  358.         ComputePoints(&Mercator,init,&grid,pout);
  359.         ComputePoints(&Mercator,&frameparam,&frame,pout);
  360.         ComputePoints(&Mercator,init,pin,pout);
  361.         break;
  362.     case MOLLWEIDE:  
  363.         ComputePoints(&Mollweide,init,&grid,pout);
  364.         ComputePoints(&Mollweide,&frameparam,&frame,pout);
  365.         ComputePoints(&Mollweide,init,pin,pout);
  366.         break;
  367.     case ORTHO: 
  368.         ComputePoints(&Ortho,init,&grid,pout);
  369.         drawCircleFrame(init->radius,pout);
  370.         ComputePoints(&Ortho,init,pin,pout); 
  371.         break;
  372.     case SINUSOIDAL: 
  373.         ComputePoints(&Sinusoidal,init,&grid,pout);
  374.         ComputePoints(&Sinusoidal,&frameparam,&frame,pout);
  375.         ComputePoints(&Sinusoidal,init,pin, pout);
  376.         break;
  377.     case STEREO: 
  378.         ComputePoints(&Stereo,init,&grid,pout);
  379.         drawCircleFrame(init->radius*1.2, pout);
  380.         ComputePoints(&Stereo,init,pin,pout);
  381.         break;
  382.     case NONE:
  383.         ComputePoints(&None,init,&grid,pout);
  384.         ComputePoints(&None,init,pin,pout);
  385.         break;
  386.     default: break;
  387.     }
  388. return 1;
  389. }
  390.  
  391.  
  392.